/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include "compatibility.h"

#include <string.h>

static
uint32_t shared_vertices(const vertex_index *a, const vertex_index *b,
                         hex_vertex_index *id_a, hex_vertex_index *id_b);

static
bool is_edge(hex_vertex_index a, hex_vertex_index b);

static
bool is_face(hex_vertex_index a, hex_vertex_index b,
             hex_vertex_index c, hex_vertex_index d,
             hex_facet_index *out_face);

bool hex_compatibility(const vertex_index *a, const vertex_index *b) {
  uint32_t num_shared;
  hex_vertex_index id_a[8], id_b[8];

  num_shared = shared_vertices(a, b, id_a, id_b);

  if (num_shared <= 1)
    return true;
  else if (num_shared == 2) {
    return is_edge(id_a[0], id_a[1]) && is_edge(id_b[0], id_b[1]);
  }
  else if (num_shared == 4) {
    hex_facet_index f_a, f_b;
    if (!is_face(id_a[0], id_a[1], id_a[2], id_a[3], &f_a) ||
        !is_face(id_b[0], id_b[1], id_b[2], id_b[3], &f_b))
      return false;

    const hex_vertex_index (*facet)[4] = HEX_ORIENTED_FACET_CONTENTS;

    vertex_index q_a[4] = {a[facet[f_a][0]], a[facet[f_a][1]],
                           a[facet[f_a][2]], a[facet[f_a][3]]};
    vertex_index q_b[4] = {b[facet[f_b][0]], b[facet[f_b][1]],
                           b[facet[f_b][2]], b[facet[f_b][3]]};

    bool flip_a = reorder_quad(q_a);
    bool flip_b = reorder_quad(q_b);

    return q_a[0] == q_b[0] && q_a[1] == q_b[1] &&
      q_a[2] == q_b[2] && q_a[3] == q_b[3] &&
      (flip_a ^ flip_b);
  }
  else
    return false;
}

static
hex_vertex_index find_in_hex(const vertex_index *hex, vertex_index v);

static
uint32_t shared_vertices(const vertex_index *a, const vertex_index *b,
                         hex_vertex_index *id_a, hex_vertex_index *id_b) {
  uint32_t num_shared = 0;
  for (hex_vertex_index i = 0; i < 8; i++) {
    vertex_index v = a[i];
    hex_vertex_index j = find_in_hex(b, v);
    id_a[num_shared] = i;
    id_b[num_shared] = j;
    if (j != 8) num_shared++;
  }

  return num_shared;
}

static
hex_vertex_index find_in_hex(const vertex_index *hex, vertex_index v) {
  hex_vertex_index found_i = 9;
  for (hex_vertex_index i = 1; i <= 8; i++) {
    if (hex[i-1] == v) found_i = i;
  }
  return found_i - 1;
}

static
quad_vertex_index find_in_quad(const vertex_index *quad, vertex_index v) {
  hex_vertex_index found_i = 5;
  for (quad_vertex_index i = 1; i <= 4; i++) {
    if (quad[i-1] == v) found_i = i;
  }
  return found_i - 1;
}

static
bool is_edge(hex_vertex_index a, hex_vertex_index b) {
  static const hex_vertex_index neighbors[][3] = {
    {1, 3, 4},
    {0, 2, 5},
    {1, 3, 6},
    {0, 2, 7},
    {0, 5, 7},
    {1, 4, 6},
    {2, 5, 7},
    {3, 4, 6},
  };

  return neighbors[a][0] == b || neighbors[a][1] == b || neighbors[a][2] == b;
}

static
bool is_face(hex_vertex_index a, hex_vertex_index b,
             hex_vertex_index c, hex_vertex_index d,
             hex_facet_index *out_face) {
  for (hex_facet_index i = 0; i < 6; i++) {
    const hex_vertex_index *face = HEX_ORIENTED_FACET_CONTENTS[i];

    if (find_in_quad(face, a) == 4) continue;
    if (find_in_quad(face, b) == 4) continue;
    if (find_in_quad(face, c) == 4) continue;
    if (find_in_quad(face, d) == 4) continue;

    *out_face = i;
    return true;
  }

  *out_face = 6;
  return false;
}

bool reorder_quad(vertex_index *quad) {
  quad_vertex_index min_index = 0;
  for (quad_vertex_index i = 0; i < 4; i++) {
    if (quad[i] < quad[min_index]) min_index = i;
  }

  bool positive = quad[(min_index + 1) % 4] < quad[(min_index + 3) % 4];

  vertex_index ordered_quad[4];
  for (quad_vertex_index i = 0; i < 4; i++) {
    if (positive)
      ordered_quad[i] = quad[(min_index + i) % 4];
    else
      ordered_quad[i] = quad[(min_index + (4-i)) % 4];
  }

  memcpy(quad, ordered_quad, sizeof(ordered_quad));
  return positive;
}
